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NUMERICAL  MODELLING  OF  SEISMIC  INTERFACE  WAVES 

by 

Henrik  Schmidt 


ABSTRACT 

s' 

/ 

/ 

SS  A  new  numerical  model;  based  on  the  linear  theory  of  elasticity  and  Marsh's 
approximation  of  the  Hankel  transform,,  has  been  developed  and  used  to  model 
the  excitation  and  propagation  of  interface  waves  in  shallow-water  sea- 
beds.  Using  realistic  material  parameters,  the  geometry  and  frequency- 
dependent  propagation  characteristics  are  analyzed  and  synthetic 
seismograms  are  produced.  Basic  properties  of  interface  waves  are 
discussed,  and  methods  for  direct  evaluation  of  sea-bed  properties  from 
experimental  results  are  proposed.  -^7 


INTRODUCTION 


The  importance  of  the  shear  properties  of  the  sea-bed  for  acoustic  wave 
propagation  in  shallow  water  is  well  established  <1>.  Unfortunately  the 
shear  parameters  are  very  difficult  to  isolate  experimentally.  They  are, 
however,  indirectly  present  through  the  properties  of  the  measurable 
seismic  interface  waves,  of  which  SACLANTCEN  has  made  considerable  experi¬ 
mental  investigation  in  recent  years  <2,3>. 

The  mathematical  model  used  in  interpreting  the  experimental  data  has 
usually  been  based  on  the  theory  of  elasticity.  Since  closed-form  solu¬ 
tions  cannot  be  obtained,  numerical  models  based  on  the  Thomson-Haskell 
technique  have  been  used.  The  shear  parameters  are  determined  indirectly 
by  fitting  the  model  predictions  to  the  experimental  data  by  parametric 
studies.  Due  to  the  comprehensive  computations  needed  at  each  frequency 
considered,  only  a  few  attempts  have  been  made  to  determine  the  shear  pro¬ 
perties  in  this  way,  <4>  for  example.  The  present  paper  describes  a  new, 
fast,  numerical  model  that  has  been  used  to  clarify  some  of  the  basic  prin¬ 
ciples  Involved  in  the  propagation  of  interface  waves  in  a  simplified 
shallow-water  environment.  The  influence  of  sediment  thickness,  frequency, 
and  material  properties  is  analyzed,  and  some  guide  lines  are  given  for 
evaluating  the  shear  properties  directly  from  the  experimental  data. 
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1  THE  MATHEMATICAL  MODEL 

The  mathematical  model  is  based  on  the  assumption  that  the  water  column  and 
the  bottom  consist  of  a  series  of  range -independent  layers.  All  materials 
are  considered  to  be  homogeneous  and  isotropic  elastic  continua  with  Lam§ 
constants  \n  and  yn  and  density  pn.  The  subscript  refers  to  layer  number 
n.  The  damping  mechanisms  are  assumed  to  be  linear  viscoelastic. 

A  cylindrical  coordinate  system  {r,e,z}  is  introduced,  with  the  z-axis  going 
through  the  source  and  being  positive  downwards  (Fig.  1).  The  representation  of 
the  cylindrical  displacement  components  {u,v,w}  in  terms  of  scalar  poten¬ 
tials  and  the  subsequent  expression  of  these  as  Hankel  transforms  closely 
follows  the  presentation  given  by  Schmidt  and  Krenk  <5>;  hence,  only  an 
outline  will  be  given  here.  If  body  forces  are  neglected,  the  displacement 
equation  of  motion  will  be  satisfied  if  the  displacement  components  in 
layer  n  are  expressed  in  terms  of  three  scalar  potentials  {*n»l'n»An}  as 


vl  n 


n 


»*n  1  1  3*n  +  »2An 
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where  the  potentials  satisfy  the  wave  equations 
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(Eq.  1) 


(Eq.  2) 


(V2 


in  which  C 
respective 


(Tn»An)  =  0.  (Eq.  3) 

cTn 

and  C-r  are  the  velocities  of  the  compressional  and  shear  waves 

y: 

^n  2un  (Eq.  4) 

Pn 


» 


Pn 

Pn 


(Eq.  5) 


In  the  present  case  the  field  is  asymmetric  due  to  the  positioning  of  the 
source  on  the  axis,  and  the  angular  displacement  v  vanishes  everywhere.  It 
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FIG.  1  SIMPLIFIED  SEA-BED  MODEL 
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Is  then  clear  from  Eq.  1  that  the  potential  Yn  must  be  constant  and  can  be 
excluded. 

In  the  following,  only  vibrations  with  angular  frequency  u  will  be 
considered;  displacements,  stresses  and  potentials  can  then  be  expressed  in 
complex  form  with  the  common  factor  exp(iut).  This  factor  will  not  be 
included  in  the  following.  The  viscoelastic  damping  can  now  be  accounted 
for  by  allowing  the  LamS  constants  to  be  complex.  After  use  of  the  Hankel 
transform  on  the  wave  equations,  the  following  integral  representations 
are  obtained  for  the  potentials: 


*n(r,z)  -/  {A-(s)e'Z°n(s)  +  A+(s)eZan(s)}  J0(rs)s  ds  (Eq.  6) 

o  n  n 

An(r,z)  =/  {B-(s)e'ZPn(s)  +  B+(s)eZB"(s)}  J0(rs)s  ds,  (Eq.  7) 

o  n  n 

where 

Jm  is  the  Bessel  function  of  the  first  kind  and  order  m. 

A",  A+,  Bjj  and  B+  are  arbitrary  functions  in  the  horizontal  wavenumber  s. 


°n(s),  Bn(s)  are  defined  as 

s  >  Re{ hn) 

«n(s)  *  J  - - 

i  /(hj-s*). 

s  <  Re{hn) 

(Eq.  8) 

(  As2-kjj, 

s  >  Re{kn) 

Pn(  ]  I  i  /(kj-s2). 

s  <  Re{kn) 

(Eq.  9) 

The  wavenumbers  hn  and  kn 
tively,  are  defined  by 

h2  n  (_“_)2  .  “2p» 

n  cLn  Xn+2yn 

kl m  -  — 

n  cTn  »*n 


for  the  compress ional 


and  shear  waves,  respec 

(Eq.  10) 

(Eq.  11) 


If  Eqs.  6  and  7  are  Inserted  into  Eq.  1,  the  following  expressions  are 
obtained  for  the  particle  displacements: 
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«(r,z)  |  A-  «■“"  ♦  A+  e“n 

o  n  n 

+  sB”  e  ZBn  +  s  B+  eZ8n}  s  J0(rs)ds,  (Eq.  12) 

n  n 

u(r,z)  I  n  “  M*s  A"  e_Z°n  ’  s  eZ°n 

+  3n  B-  e”Z3n  -  Bn  B+  eZB"}  s  Ji(rs)ds.  (Eq.  13) 

n  n 

The  stress  components  involved  in  the  boundary  conditions  follow,  by 
Hooke's  law: 

Ozzlr.*)  |  „  ■  u„  1  ((2s2-k2)(A-  e"zan  ♦  AJ  e“") 

0 

+  2sBn(-B-  e’z6n  +  B+  ez8n)}  s  J0(rs)ds,  (Eq.  14) 
n  n 

®rz(r»z)  |  n  =  Wn  /  {2san(Aj  e  -  A+  e  n) 

-(2sz-k2)(B_  e‘ZBn  +  B+  eZBn)}  s  J,(rs)  ds. 

n  n  n  (Eq.  15) 

In  the  case  of  a  fluid  layer  the  shear  stiffness  un  vanishes,  and  only  the 
potential  *n  is  involved.  The  displacements  follow  directly  from  Eqs.  12 
and  13  by  setting  B"nand  B+nto  zero.  The  shear  stress  is  identically  zero, 

whereas  Eq.  14  has  to  be  replaced  by 

•  -zan 

«zz(r,z)  |  n  ■  -M2  /  1A‘  e  +  AJ  e6Zan}  s  J0(rs)ds.  (Eq.  16) 


The  source  is  assumed  to  be  in  layer  number  m  at  depth  zs.  In  the  absence 
of  boundaries,  the  field  produced  in  layer  m  would  be  <6>: 


♦s(r.z) 


jSg, 

4w 


oo  -  |  Z-ZS  |  am 

/  - - s  J0(rs)ds 

0  «m 


As(r,z)  s  0  , 


(Eq.  17) 
(Eq.  18) 


where  Is  the  source  strength.  If  Eq.  17  Is  Inserted  into  Eq.  1, 
expressions  similar  to  Eqs.  12  and  13  are  obtained  for  the  displacements. 
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and  again  Hooke's  law  yields  expressions  like  Eqs.  14  and  IS  for  the 
stresses  Involved  In  the  boundary  conditions. 

For  each  value  of  the  range,  r,  the  boundary  conditions  must  be  satisfied. 
In  the  upper  and  lower  half-spaces  the  arbitrary  functions,  with 
superscript  •  and  +  respectively,  must  vanish  due  to  the  radiation  con¬ 
dition.  At  each  interface,  w  and  azz  must  be  continuous,  and  at  all 
solld/llquld  Interfaces  the  shear  stresses  must  vanish.  At  solld/solld 
Interfaces  w,  u,  oZ2  and  orz  must  be  continuous.  This  yields  a  linear 
system  of  equations  in  the  arbitrary  functions,  to  be  satisfied  at  each 
horizontal  wavenuriber: 

Cij(s)  •  Aj(s)  *  R-j (s)  .  (Eq.  19) 

The  vector  Aj(s)  contains  all  the  non-vanishing  arbitrary  functions,  C-j j (s) 
is  the  coefficient  matrix,  and  R<j (s)  contains  the  contributions  from  the 
source.  When  the  arbitrary  functions  are  found,  the  field  parameters  at 
any  depth  and  range  can  be  obtained  from  the  Hankel  transforms  (Eqs.  12  to 
16)  plus  the  source  contributions  (if  the  source  and  receiver  are  in  the 
same  layer). 

An  analytical  solution  of  Eq.  19  is  of  course  possible,  leading  to  closed- 
form  expressions  for  the  arbitrary  functions;  but  for  more  than  a  few 
layers  this  procedure  would  be  inconvenient.  Further,  the  Hankel  trans¬ 
forms  do  not  lead  to  closed-form  solutions,  but  need  numerical  evaluation. 
Thus  the  most  general  way  to  proceed  is  to  create  a  numerical  model  based 
directly  on  the  system  of  equations  (Eq.  19).  Such  a  model  is  described  in 
the  next  chapter. 


2  THE  NUMERICAL  MODEL 

The  numerical  evaluation  of  the  Hankel  transform  necessitates  a  truncation 
and  a  discretization  in  the  horizontal  wavenumber  s.  As  can  be  observed 
from  Eqs.  17  and  8,  the  source  terms  decay  exponentially  for  s  going  towards 
infinity.  As  the  source  terms  form  the  right  hand  side  of  Eq.  19  the 
arbitrary  functions  will  behave  in  the  same  way.  It  is  therefore  possible 
to  truncate  the  Integration  interval  in  accordance  with  any  accuracy 
demands.  The  fast-field  technique  Introduced  by  Marsh  <7>  can  then  be 
used  to  evaluate  the  Hankel  transforms. 

The  Bessel  functions  are  expressed  in  terms  of  Hankel  functions 

Jm(rs)  *  *  Hil)<rs>  +  Hi2)(rs)  (Eq*  W) 

and  each  Integral  Is  split  into  two.  As  only  outgoing  waves  are  considered, 
the  Integrals  involving  H(l)  (rs)  are  neglected,  and  Hjr)  (rs)  is  replaced  by 
Its  asymptotic  form  m 


(rs) 


'V. 

rs-*» 


2  i 

<*»> 


-1[rs-(m+J)  £] 
e 


(Eq.  2D 
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The  integration  over  the  truncated  interval  can  now  be  performed  by  means 
of  the  fast  fourier  transform,  and  the  actual  field  parameter  is  found  at  a 
number  of  ranges  equal  to  the  number  of  discrete  wavenumbers  considered. 

DiNapoli  and  Deavenport  <8>  have  compared  Marsh's  method  to  the  technique 
introduced  by  Tsang,  Brown,  Kiang  and  Simeons  <9>,  which  does  not  use  the 
assymptotic  form  of  the  Hankel  function  and  found  significant  differences 
only  for  very  short  ranges. 

The  kernels  in  the  Hankel  transforms  are  now  needed  only  for  a  limited 
number  of  discrete  values  of  s.  Kutschale  <10>  and  DiNapoli  and  Deavenport 
<8>  used  a  Greens-function  approach  based  on  the  Thomson-Haskel 1  matrix 
method.  However,  their  approach  allows  for  only  one  source/ receiver  com¬ 
bination  at  a  time.  Here  we  solve  Eq.  19  directly.  The  receiver  depth  is 
not  involved  in  the  equations  and  several  receiver  depths  can  therefore  be 
handled  with  one  solution.  Further,  the  coefficient  matrix  does  not 
include  the  source  contribution,  thus  yielding  the  possibility  of  having 
several  sources  on  the  axis  by  simply  adding  their  contributions  to  form 
the  righthand  side  of  Eq.  19. 

The  solution  of  Eq.  19  is  the  most  critical  point  with  respect  to  com¬ 
putation  time,  but  if  only  the  non-vanishing  arbitrary  functions  are 
included  and  arranged  properly,  the  coefficient  matrix  will  be  compact  and 
of  band  form.  The  equations  are  then  solved  very  efficiently  by  means  of 
gaussian  elimination  with  partial  pivoting.  In  all  cases  treated,  the 
total  calculation  time  with  this  model  was  between  10  and  30%  of  what  was 
used  by  Kutschale's  model,  with  the  biggest  gain  in  cases  with  many  layers. 
In  cases  with  more  than  one  source  and  receiver  the  gain  is  of  course  much 
more  pronounced. 


3  MODELLING  OF  SEISMIC  WAVES 

The  numerical  model  has  been  used  to  analyze  the  propagation  of  low- 
frequency  seismic  waves  in  shallow-water  sea-beds.  In  order  not  to  obscure 
the  basic  principles,  a  simple  two-layered  model  was  chosen  for  the  sea¬ 
bed,  Fig.  1.  Below  the  water  column  of  depth  d#,  a  single  sediment  layer 
of  thickness  dc  covers  a  half-space  of  rock  or  rock-like  material.  The 
bottom  materials  used  in  the  examples  and  their  assumed  properties  are 
listed  in  Table  1.  The  complex  Lam§  constants  are  not  given  explicitly  in 
the  table.  Instead,  the  compressional  and  shear  velocities  are  shown, 
together  with  their  respective  dampings  in  decibels  per  wavelength. 


( 
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TABLE  1 

MATERIAL  PROPERTIES 


Table  2  summarizes  the  test  cases  considered.  As  can  be  observed,  the 
water  depth  is  held  constant  in  all  cases.  Soft  and  hard  sediment  layers 
of  50  m  thickness  are  combined  with  two  different  sub-bottom  materials  in 
order  to  assess  the  effect  of  the  material's  properties.  For  one  of  the 
combinations  the  thickness  of  the  sediment  layer  is  varied.  In  the  last 
test  case  the  effect  of  a  shear-speed  gradient  in  the  sediment  layer  is 
analyzed  by  subdividing  the  50-m  layer  into  ten  homogeneous  5-m  layers. 


TABLE  2 
TEST  CASES 


Test  Case 

Water 

Depth 

(m) 

Sediment 

Thickness 

On) 

Sediment 

Material 

shear-Speed 

Gradient 

(m/s/m) 

Sub-bottom 

Material 

1 

100 

50 

Silt 

. 

Limestone 

2 

100 

30 

Silt 

- 

Limestone 

3 

100 

5 

Silt 

- 

Limestone 

4 

100 

50 

Sand 

• 

Limestone 

5 

100 

50 

Silt 

• 

Basalt 

6 

100 

50 

Sand 

- 

Basalt 

7 

100 

50 

Silt 

4 

Limestone 

In  all  cases  a  point  source  is  placed  in  the  middle  of  the  water  column  at 
50  m  depth,  and  the  vertical  particle  velocity  at  the  top  of  the  sediment 
layer  is  calculated  using  a  pressure  amplitude  of  IPa  at  a  distance  of  1  m 
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from  the  source.  Only  frequencies  below  and  around  the  cut-off  frequency 
for  the  first  water  mode  are  considered,  which  in  the  present  case  means 
frequencies  below  10  Hz. 

In  the  first  test  case  a  50-m  silt  sediment  layer  covers  a  limestone  sub¬ 
bottom.  Figure  2  shows  the  modulus  of  the  integrand  in  the  Hankel  trans¬ 
form  of  the  vertical  velocity  at  the  top  of  the  sediment  layer  for  a 
frequency  of  1.5  Hz.  No  water  modes  are  present  at  this  frequency,  but  two 
mode-like  peaks  can  be  observed  at  wavenumbers  corresponding  to  phase  velo¬ 
cities  of  300  m/s  (peak  1)  and  840  m/s  (peak  2).  These  peaks  correspond  to 
interface  waves,  and  will  be  denoted  the  first  and  second  interface  mode, 
respectively.  The  first  interface  mode  is  best  excited  (highest 
amplitude),  but  it  also  has  the  highest  damping  (widest  peak). 

Figure  3  shows  the  corresponding  vertical  particle  velocities  at  ranges  up 
to  50  km.  Due  to  the  high  damping  of  the  first  interface  mode,  its  contri¬ 
bution  is  significant  only  for  ranges  shorter  than  2  km,  beyond  which  the 
second  interface  mode  becomes  dominant. 

If  the  source  is  not  of  stationary  type,  but  transient,  the  different  velo¬ 
cities  of  the  two  interface  modes  will  yield  different  arrival  times,  and 
the  presence  of  the  first  interface  mode  will  be  measurable  also  at  greater 
ranges.  The  phase  and  group  velocities  of  the  two  interface  modes  have 
been  calculated  in  the  frequency  band  of  interest,  0.1  to  10  Hz,  with  a 
resolution  of  0.1  Hz.  The  results  are  shown  in  Fig.  4  together  with  an 
excitation  measure,  which  somewhat  arbitrarily  has  been  chosen  to  be  the 
particle  velocity  at  10  km  range. 

Since  most  of  the  general  features  of  Fig.  4  exist  in  all  test  cases,  a 
detailed  interpretation  will  be  given  at  this  point.  There  is  only  one 
interface  mode  at  frequencies  below  1  Hz.  It  is  slightly  dispersive,  with 
phase  and  group  velocities  approaching  those  of  a  Rayleigh  wave  on  a 
limestone  half-space.  At  1  Hz  a  very  sharp  transition  zone  appears,  where 
the  velocities  drop  dramatically  to  values  approaching  those  of  a  Scholte 
wave  at  a  water/silt  interface.  At  1.8  Hz  the  group  velocity  reaches  a 
minimum  of  100  m/s,  i.e.  half  the  shear  speed  in  silt.  The  second  inter¬ 
face  mode  has  a  sharp  cut-off  at  the  transition  frequency,  and  after  a 
distinct  minimum  of  the  group  velocity  it  appears  as  a  logical  continuation 
of  the  low-frequency  part  of  the  first  interface  mode.  Above  2  Hz  the 
excitation  (solid  line)  of  the  second  interface  mode  decreases  due  to  the 
increased  distance  in  terms  of  wavelengths  of  the  source  from  the 
silt/limestone  interface  along  which  the  second  interface  mode  propagates. 
The  sharp  peak  on  the  excitation  curve  around  4.5  Hz  is  the  first  propaga¬ 
tion  water  mode. 

The  existence  of  mode  transition  zones  is  well  known  from  the  theory  of 
vibration  of  elastic  plates,  Mindlin  <11>,  where  they  appear  near  the 
thickness-shear  frequencies.  These  are  the  frequencies  at  which  an  infi¬ 
nite  elastic  plate  can  perform  free  shear  vibrations  with  vanishing  ver¬ 
tical  displacements.  Now  consider  the  silt  layer  as  an  infinite  elastic 
plate  in  welded  contact  with  an  infinite  rigid  half-space.  The  first 
thickness-shear  frequency  would  then  correspond  to  that  of  a  free  silt 
plate  of  the  double  thickness,  <11>: 


VERTICAL  PARTICLE  VELOCITY  (dB//1m/s) 
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FIG.  2  INTEGRAND  FOR  VERTICAL  PARTICLE  VELOCITY  AT  1.5  Hz 
TEST  CASE  1:  50  a  SILT  ON  LIMESTONE 


FIG.  3  VERTICAL  PARTICLE  VELOCITIES  AT  1.5  Hz 
TEST  CASE  1:  SO  m  SILT  ON  LIMESTONE 
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fTS  *  Ta  .  (Eq.  22) 

4ds 

where  Cj  is  the  shear  velocity  and  ds  is  the  thickness  of  the  layer.  If 
the  parameters  for  silt  are  applied  to  Eq.  1,  we  obtain  fyj  equal  to  1  Hz, 
which  is  very  close  to  the  observed  transition  frequency  in  Fig.  4. 

Since  the  transition  frequency  can  be  observed  in  experimental  results, 
<3>,  its  correlation  with  the  thickness-shear  frequency  could  yield  a 
direct  method  of  determining  the  shear-wave  velocity  in  a  single  sediment 
layer  overlying  a  rigid  half-space. 

To  summarize  the  general  propagation  characteristics  for  test  case  1 
(Fig.  4),  the  computed  particle  velocity  at  10  km  range  is  entirely  asso¬ 
ciated  with  the  first  interface  mode  below  the  transition  frequency  (1  Hz). 
This  interface  mode  is  strongly  related  to  a  Rayleigh  wave  on  a  limestone 
half-space  below  1  Hz,  while  it  becomes  an  interface  wave  connected  with 
the  water/silt  interface  at  frequencies  above  1  Hz.  In  this  frequency 
regime  the  second  interface  mode  appears,  and  it  is  mainly  related  to  the 
silt/limestone  interface.  Finally,  the  first  water-borne  mode  appears  at 
around  4.5  Hz. 

In  test  case  2  (Fig.  5)  the  thickness  of  the  silt  layer  has  been  reduced  to 
30  m.  The  general  behaviour  of  the  phase  and  group  velocities  is  the  same 
as  in  the  former  example,  but  the  transition  frequency  is  now  at  approxima¬ 
tely  1.6  Hz,  which  is  in  good  agreement  with  the  computed  thickness-shear 
frequency  of  1.67  Hz.  If  the  thickness  of  the  silt  layer  is  further 
reduced  to  5  m,  the  thickness-shear  frequency  becomes  10  Hz,  and  no  tran¬ 
sition  zone  is  present  in  the  frequency  range  of  interest. 

In  the  fourth  test  case,  the  silt  sediment  is  replaced  by  a  sand  layer  of 
50  m  thickness.  The  transition  zone  is  less  sharp  in  this  case  (Fig.  6), 
probably  due  to  the  smaller  difference  between  the  stiffness  of  the  sedi¬ 
ment  and  the  sub-bottom.  In  this  case,  the  thickness-shear  frequency  is 
3  Hz  and  the  minimum  group  velocity  is  400  m/s,  or  67%  of  the  shear  speed. 
At  the  low-frequency  end  the  velocities  of  the  first  interface  mode  again 
approach  the  velocity  of  a  Rayleigh  wave  on  the  surface  of  a  limestone 
half-space. 

Silt  and  sand  sediment  layers  (50  m)  on  basalt  were  also  studied  in  order 
to  analyze  the  effect  of  different  sub-bottom  materials.  In  the  case  of 
silt  on  basalt  (Fig.  7)  the  low-frequency  part  of  the  first  interface  mode 
is  excited  only  slightly,  but  the  transition  zone  can  again  be  observed 
near  the  thickness-shear  frequency  of  1  Hz.  Above  1  Hz  the  dispersion  cur¬ 
ves  are  very  similar  to  those  obtained  earlier  for  a  sub-bottom  of 
limestone  (Fig.  4).  In  the  case  of  sand  on  basalt  (Fig.  8),  the  effect  of 
the  sub-bottom  material  is  more  pronounced.  Compared  with  the  case  for 
limestone  (Fig.  6),  the  basalt  sub-bottom  yields  a  much  sharper  transition 
zone  near  the  thickness-shear  frequency  of  the  sand  layer  (3  Hz).  The 
group  velocity  has  a  minimum  of  one-half  of  the  shear  speed  in  sand  and  the 
Rayleigh-wave  behaviour  at  low  frequencies  can  again  be  observed. 

To  Investigate  the  effect  of  a  shear  gradient  in  the  sediment  layer,  test 
case  1  (Fig.  4)  was  modified  to  Include  a  gradient  of  4  m/s/m  in  the  silt 
layer.  The  result  Is  given  In  Fig.  9.  The  main  effect  Is  seen  to  be  a 
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FIGS.  4  to  9  EXCITATION  AND  DISPERSION  CURVES 
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TEST  CASE  1:  50  m  SILT  ON  LIMESTONE 
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TEST  CASE  2:  30  m  SILT  ON  LIMBSTON 


c o  1000 


5 

g,eoo 

5^  400 
1 

»  200 


0.1  1  to 

FREQUENCY  (Hz) 

FIG.  6 

TEST  CASE  4:  50  m  SAND  ON  LIMESTONE 


-200  -S 

c/3  1000 

? 

-220  ^ 

1 

§  800 
-f 

£ 

-240  £ 

%^600 

O 

g  « 

-260  ^ 

P4oo 

V 

-280  ^ 

i 

H  200 

5 

-300  *■ 

£  0 

to  2000 
S  1000 


FREQUENCY  (Hz) 

FIG.  8 

TEST  CASE  6:  SO  m  SAND  ON  BASALT 
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TEST  CASE  5:  50  m  SILT  ON  BASALT 
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TEST  CASE  7:  50  n  SILT  MITE  SHEAR 

SPEED  GRADIENT  ON  LIMESTONE 
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change  in  the  high-frequency  limit  of  the  dispersion  curves  for  the  first 
interface  mode.  Since  the  average  shear-wave  velocity  is  the  same  as  in 
test  case  1,  the  transition  frequency  is  not  affected.  The  effect  on  the 
high-frequency  limit  is  not  very  surprising.  The  shorter  the  wavelength, 
the  more  dominant  is  the  upper  part  of  the  sediment,  yielding  lower  phase 
and  group  velocities.  The  almost  non-dispersive  nature  of  the  first  inter¬ 
face  mode  in  this  case  is  often  seen  in  experimental  data,  <3>. 


4  TIME  SERIES 

Although  no  attempt  has  been  made  to  model  experimental  results,  synthetic 
seismograms  have  been  produced  for  two  of  the  test  cases  in  order  to 
illustrate  the  time-domain  effect  of  the  features  described  above.  The 
source  is  assumed  to  be  half  a  sine  wave  of  1.5  Hz  sent  through  an  ideal 
0.5  to  3.2  Hz  band-pass  filter.  This  frequency  range  has  been  chosen 
because  it  contains  frequencies  on  both  sides  of  the  transition  frequency 
for  a  silt  sediment  layer  of  50  m  thickness.  The  transfer  functions  were 
calculated  to  a  resolution  of  0.01  Hz  and  multiplied  by  the  spectrum  of  the 
source.  The  time  series  were  then  created  by  means  of  the  fast  fourier 
transform  at  ranges  of  1,  2,  3,  4  and  5  km. 

Figure  10  shows  the  result  for  test  case  1,  i.e.  a  50  m  silt  layer  on  a 
limestone  half-space.  The  first  weak  arrival  corresponds  to  a 

compressional  wave  in  the  limestone,  whereas  the  first  significant  arrival 
corresponds  to  the  Rayleigh  wave  velocity  of  the  limestone  (900  m/s).  This 
arrival  consists  of  the  low-frequency  parts  of  the  first  and  second  inter¬ 
face  modes.  A  clear  dispersion  can  be  observed  corresponding  to  the 
distinct  minimum  in  the  group  velocity  of  the  second  interface  mode  in 

Fig.  4.  The  slow,  highly  damped  wave,  corresponding  to  the  first  interface 
mode  above  the  transition  frequency,  propagates  with  group  velocities  bet¬ 
ween  100  and  180  m/s,  again  in  good  agreement  with  Fig.  4.  At  ranges 

greater  than  4  km  this  arrival  has  negligible  amplitude. 

As  described  above,  the  low-frequency  part  of  the  first  interface  mode  will 
not  be  significantly  excited  if  the  sub-bottom  is  basalt.  The  maximum 
excitation  at  10  km  range  lies  at  2.5  Hz  (Fig.  7)  and  is  due  to  a  strong 
excitation  of  the  second  mode.  These  properties  are  also  reflected  in  the 
synthetic  seismograms  (Fig.  11).  The  fastest  arrival,  apart  from  the  weak 
compressional  wave,  has  its  major  frequency  content  above  2  Hz  and 

corresponds  mainly  to  the  high-frequency  end  of  the  second  interface  mode. 
The  severe  dispersion  of  the  corresponding  arrival  in  Fig.  10  is  not  pre¬ 
sent  here,  and  the  slow  part  of  the  first  interface  mode  is  well  separated, 
travelling  at  group  velocities  between  100  and  200  m/s.  The  synthetic 
seismograms  in  Figs.  10  and  11  are,  at  least  qualitatively,  very  similar  to 
those  observed  during  experiments,  <2,3>. 
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STACKED  SYNTHETIC  SEISMOGRAMS 
TEST  CASE  1:  50  m  SILT  ON  LIMESTONE 
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CONCLUSIONS 


A  fast  numerical  model  based  on  the  linear  theory  of  elasticity  has  been 
developed  and  used  for  the  modelling  of  seismic  waves  in  a  shallow-water 
sea-bed.  It  is  shown  that  a  model  of  this  type  can  be  used  to  analyze  the 
excitation  and  propagation  of  interface  waves,  reproducing  many  of  the 
basic  features  observed  during  experiments,  including  the  presence  of  a 
characteristic  group-velocity  transition  zone. 

One  of  the  main  prospects  of  the  investigation  of  interface  waves  is  the 
possibility  of  determining  material  properties  that  are  not  directly 
measurable,  such  as  shear-wave  velocities  and  dampings.  Below  the  tran¬ 
sition  frequency,  the  propagation  of  the  first  interface  mode  is  similar  to 
that  of  a  Rayleigh  wave  on  a  pure  sub-bottom  half-space.  Its  velocity  and 
damping  therefore  yield  information  on  the  properties  of  the  sub-bottom. 
The  transition  frequency  is  closely  related  to  the  first  thickness-shear 
frequency  of  the  sediment  layer  and  is  therefore  a  direct  measure  of  the 
sediment  shear  speed.  Above  the  transition  frequency,  the  first  interface 
mode  is  closely  related  to  the  water/sediment  interface,  and  its  properties 
can  therefore  be  used  to  determine  the  properties  of  the  sediment  layer. 
The  second  interface  mode,  appearing  above  the  transition  frequency,  is 
related  to  the  sediment/sub-bottom  interface  and  therefore  contains  coupled 
information  on  both  media.  If  more  than  one  sediment  layer  is  present, 
several  interface  modes  will  appear,  and  the  interpretation  of  the  theore¬ 
tical  results  can  become  quite  complex. 
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